#!/usr/bin/perl

# ------------------------------------------------------------------------------
# --- Get FASTA from DSSP secondary structure file (ss.txt) for each PDB ID ----
# --- http://www.rcsb.org/pdb/static.do?p=download/http/index.html -------------
# --- http://www.rcsb.org/pdb/files/ss.txt.gz ----------------------------------
# --- Original by A. Jochim ----------------------------------------------------
# --- Modified and documented by C. Bergey and A. Watkins ----------------------
# ------------------------------------------------------------------------------

use strict;
use warnings;

# Be sure to run "module load perl" from the PBS script first.
use Bio::DB::Fasta;
use Bio::SeqIO;

print STDERR "Indexing FASTA file of secondary structures and sequences...\n";

# First, open latest list of new PDB IDs
my $latest_pdb_ids = shift;
open PDB_IDS, ("<$latest_pdb_ids")
	or die "Could not open file of new PDB IDs ($latest_pdb_ids)\n";
my @pdb_ids = <PDB_IDS>;

print STDERR "\tNumber of IDs in PDB ID file:\t", scalar @pdb_ids, "\n";

my $db = Bio::DB::Fasta->new('ss.txt');

my $seq_out = Bio::SeqIO->new(	-file =>   ">pdb_dssp_seq.fa",
                                -format => 'fasta');

# Get all IDs in the FASTA file of secondary structures and sequences
my @dssp_ids = $db->ids;
	
print STDERR "\tNumber of IDs in DSSP FASTA file:\t", scalar @dssp_ids, "\n";
	
print STDERR "Extracting FASTA sequences from secondary structure file...\n";
	
# For each new complex: 
foreach my $id (@pdb_ids) {

	chomp $id;

	# Look for DSSP IDs that match it.
	my @dssp_match_ids = grep /($id)/, @dssp_ids;
	
	# For each matching ID...
	foreach my $dssp_match_id (@dssp_match_ids) {
   
		# Grab the sequence object
		my $obj = $db->get_Seq_by_id($dssp_match_id);

		# Print it to file
		if ($obj->display_id =~ /sequence/) {
			$seq_out->write_seq($obj);
		}
      
	}
}

# Print out number of sequences and structures found in DSSP file
my $dssp_seq_count = `grep '>' pdb_dssp_seq.fa | wc -l`;
$dssp_seq_count =~ s/\s//g;

print STDERR "\tNumber of sequences found in DSSP file:\t$dssp_seq_count\n";
